stanEx5
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}hierarchical model
class c=1~k
b0c~N(b00,sb0)
b1c~N(b10,sb1)
yc~N(b0c+b1c*x,s)n=50
k=9
x=runif(n,0,20)
c=sample(letters[1:k],n,replace=T)
b00=rnorm(k,10,5)
b0=b00[factor(c)]
b10=rnorm(k,2,1)
b1=b10[factor(c)]
y=rnorm(n,b0+b1*x,5)
qplot(x,y,shape=c,size=I(2))+
scale_shape_manual(values=1:k)
qplot(x,y,col=c)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionestimate as no class
ex8-0.stan
y~N(b00+b10*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}mdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
fn(mdl,data)## Initial log joint probability = -456.673
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 23 -165.36 0.00977385 0.00137419 1 1 24
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -165.36
## b0 15.04
## b1 1.57
## s 16.56
## m0[1] 19.24
## m0[2] 30.59
## m0[3] 38.28
## m0[4] 25.01
## m0[5] 15.61
## m0[6] 18.88
## m0[7] 36.08
## m0[8] 38.06
## m0[9] 31.31
## m0[10] 39.85
## m0[11] 46.01
## m0[12] 33.20
## m0[13] 42.04
## m0[14] 31.87
## m0[15] 25.03
## m0[16] 25.33
## m0[17] 43.10
## m0[18] 32.24
## m0[19] 41.91
## m0[20] 35.20
## m0[21] 30.50
## m0[22] 21.31
## m0[23] 41.22
## m0[24] 45.78
## m0[25] 45.40
## m0[26] 43.43
## m0[27] 21.12
## m0[28] 18.46
## m0[29] 23.21
## m0[30] 25.39
## m0[31] 39.09
## m0[32] 40.21
## m0[33] 33.91
## m0[34] 23.15
## m0[35] 41.73
## m0[36] 24.18
## m0[37] 35.80
## m0[38] 43.58
## m0[39] 28.55
## m0[40] 31.78
## m0[41] 15.39
## m0[42] 42.26
## m0[43] 22.36
## m0[44] 25.72
## m0[45] 37.82
## m0[46] 40.62
## m0[47] 40.85
## m0[48] 26.78
## m0[49] 34.30
## m0[50] 22.82
## y0[1] 21.34
## y0[2] 56.24
## y0[3] 16.27
## y0[4] -9.63
## y0[5] -6.93
## y0[6] 38.76
## y0[7] 37.65
## y0[8] 43.24
## y0[9] 3.28
## y0[10] 41.76
## y0[11] 48.98
## y0[12] 1.68
## y0[13] 30.40
## y0[14] 42.67
## y0[15] -5.45
## y0[16] 53.63
## y0[17] 42.79
## y0[18] 49.27
## y0[19] 52.37
## y0[20] 44.08
## y0[21] 60.54
## y0[22] 30.18
## y0[23] 48.02
## y0[24] 45.54
## y0[25] 28.12
## y0[26] 63.27
## y0[27] 23.60
## y0[28] 45.35
## y0[29] 34.86
## y0[30] 42.53
## y0[31] 32.93
## y0[32] 30.24
## y0[33] 44.82
## y0[34] 16.25
## y0[35] 53.02
## y0[36] 19.29
## y0[37] 25.12
## y0[38] 31.93
## y0[39] 43.55
## y0[40] 31.07
## y0[41] 16.88
## y0[42] 44.27
## y0[43] 8.55
## y0[44] 25.93
## y0[45] 38.93
## y0[46] 44.69
## y0[47] 62.46
## y0[48] 59.17
## y0[49] 15.37
## y0[50] 49.14
## m1[1] 15.39
## m1[2] 18.79
## m1[3] 22.20
## m1[4] 25.60
## m1[5] 29.00
## m1[6] 32.40
## m1[7] 35.80
## m1[8] 39.20
## m1[9] 42.61
## m1[10] 46.01
## y1[1] -0.23
## y1[2] 22.26
## y1[3] 22.35
## y1[4] 32.95
## y1[5] 12.85
## y1[6] 69.10
## y1[7] 54.09
## y1[8] 31.67
## y1[9] 44.55
## y1[10] 49.48
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -164.17 -163.85 1.35 1.08 -166.81 -162.74 1.01 439 619
## b0 14.52 14.54 5.84 5.78 4.43 23.33 1.02 216 256
## b1 1.60 1.61 0.47 0.47 0.86 2.40 1.02 260 396
## s 17.38 17.24 1.87 1.82 14.57 20.61 1.00 1984 1300
## m0[1] 18.82 18.87 4.74 4.67 10.62 25.95 1.02 218 247
## m0[2] 30.42 30.38 2.54 2.46 26.08 34.44 1.01 403 656
## m0[3] 38.27 38.31 2.91 2.73 33.36 42.97 1.00 1402 1357
## m0[4] 24.71 24.80 3.39 3.25 18.73 29.89 1.01 239 279
## m0[5] 15.11 15.14 5.69 5.65 5.27 23.71 1.02 216 248
## m0[6] 18.45 18.53 4.83 4.77 10.11 25.74 1.02 217 248
## m0[7] 36.03 36.02 2.61 2.50 31.67 40.34 1.00 1459 1316
## m0[8] 38.05 38.07 2.87 2.68 33.17 42.76 1.00 1423 1325
## m0[9] 31.15 31.14 2.49 2.44 26.95 35.13 1.00 467 785
## m0[10] 39.88 40.01 3.19 3.00 34.52 45.13 1.00 1146 1403
## m0[11] 46.17 46.30 4.58 4.40 38.39 53.70 1.01 578 1338
## m0[12] 33.09 33.10 2.44 2.38 29.01 37.05 1.00 749 1049
## m0[13] 42.12 42.26 3.64 3.45 35.96 47.97 1.01 857 1363
## m0[14] 31.73 31.75 2.46 2.41 27.59 35.72 1.00 534 832
## m0[15] 24.74 24.82 3.38 3.24 18.76 29.92 1.01 239 288
## m0[16] 25.04 25.11 3.32 3.18 19.19 30.18 1.01 242 289
## m0[17] 43.21 43.34 3.88 3.65 36.67 49.47 1.01 752 1303
## m0[18] 32.11 32.12 2.45 2.38 27.97 36.10 1.00 580 865
## m0[19] 41.99 42.13 3.61 3.41 35.86 47.80 1.00 874 1363
## m0[20] 35.13 35.09 2.53 2.41 30.87 39.26 1.00 1312 1306
## m0[21] 30.33 30.29 2.54 2.47 25.97 34.39 1.01 396 656
## m0[22] 20.93 20.98 4.22 4.16 13.63 27.29 1.02 220 243
## m0[23] 41.28 41.39 3.46 3.29 35.43 46.91 1.00 964 1375
## m0[24] 45.94 46.06 4.52 4.34 38.24 53.32 1.01 588 1338
## m0[25] 45.56 45.68 4.43 4.26 38.03 52.77 1.01 605 1303
## m0[26] 43.54 43.67 3.96 3.73 36.90 49.96 1.01 725 1321
## m0[27] 20.74 20.78 4.27 4.21 13.37 27.15 1.02 220 243
## m0[28] 18.02 18.13 4.94 4.91 9.51 25.46 1.02 217 248
## m0[29] 22.88 22.98 3.78 3.68 16.32 28.57 1.02 227 257
## m0[30] 25.10 25.18 3.31 3.17 19.26 30.22 1.01 243 299
## m0[31] 39.11 39.18 3.05 2.83 34.04 44.04 1.00 1262 1357
## m0[32] 40.25 40.35 3.26 3.10 34.75 45.67 1.00 1092 1404
## m0[33] 33.81 33.80 2.46 2.38 29.69 37.77 1.00 916 1252
## m0[34] 22.82 22.91 3.79 3.70 16.25 28.52 1.02 227 257
## m0[35] 41.80 41.94 3.57 3.37 35.76 47.59 1.00 894 1363
## m0[36] 23.87 23.96 3.56 3.42 17.58 29.23 1.02 232 275
## m0[37] 35.74 35.70 2.58 2.48 31.41 39.98 1.00 1437 1347
## m0[38] 43.69 43.81 3.99 3.78 36.98 50.17 1.01 714 1321
## m0[39] 28.33 28.35 2.76 2.68 23.58 32.67 1.01 302 458
## m0[40] 31.64 31.66 2.46 2.41 27.48 35.63 1.00 522 827
## m0[41] 14.88 14.91 5.75 5.70 4.93 23.57 1.02 216 256
## m0[42] 42.35 42.50 3.69 3.48 36.12 48.27 1.01 834 1325
## m0[43] 22.01 22.09 3.97 3.88 15.21 27.98 1.02 224 245
## m0[44] 25.44 25.51 3.24 3.13 19.78 30.48 1.01 246 301
## m0[45] 37.80 37.81 2.84 2.64 32.96 42.42 1.00 1443 1326
## m0[46] 40.67 40.77 3.34 3.19 35.04 46.18 1.00 1037 1402
## m0[47] 40.90 41.00 3.39 3.23 35.19 46.45 1.00 1008 1401
## m0[48] 26.52 26.58 3.05 2.91 21.30 31.32 1.01 261 320
## m0[49] 34.21 34.20 2.47 2.41 30.10 38.20 1.00 1023 1295
## m0[50] 22.47 22.56 3.87 3.78 15.83 28.30 1.02 225 257
## y0[1] 19.28 19.53 18.18 17.46 -10.89 48.86 1.00 1528 1827
## y0[2] 30.58 31.07 17.53 17.53 0.40 58.64 1.00 1891 1957
## y0[3] 38.38 38.45 17.66 17.84 8.96 66.76 1.00 1851 1913
## y0[4] 24.63 24.80 17.78 17.50 -5.25 52.96 1.00 1793 2087
## y0[5] 15.20 15.50 18.09 17.53 -15.41 44.97 1.00 1297 1754
## y0[6] 18.21 18.47 18.41 17.62 -12.14 47.98 1.00 1573 1357
## y0[7] 36.38 35.88 17.32 16.86 8.85 65.82 1.00 1681 1895
## y0[8] 37.63 37.23 17.38 17.22 9.93 66.49 1.00 1879 1869
## y0[9] 31.34 31.12 17.93 17.53 2.66 60.82 1.00 1833 1933
## y0[10] 39.58 39.49 18.03 17.82 10.26 69.25 1.00 1977 1833
## y0[11] 46.69 46.58 18.16 17.87 15.53 75.59 1.00 1481 1368
## y0[12] 33.40 33.45 17.94 17.88 4.41 62.14 1.00 1721 1973
## y0[13] 41.98 41.97 17.53 17.23 12.71 71.17 1.00 2064 1970
## y0[14] 31.77 31.71 17.63 17.42 3.45 60.27 1.00 2199 1985
## y0[15] 25.20 25.68 17.63 17.21 -4.45 53.96 1.00 1979 1701
## y0[16] 24.97 25.30 17.57 17.19 -5.30 53.80 1.00 1220 1950
## y0[17] 43.53 43.43 18.13 17.78 13.37 72.74 1.00 1853 1747
## y0[18] 32.37 32.25 17.40 16.85 3.99 60.75 1.00 2067 1953
## y0[19] 41.41 41.05 17.68 17.49 13.66 70.50 1.00 1809 1834
## y0[20] 35.32 34.97 17.52 16.71 6.83 63.94 1.00 2012 1895
## y0[21] 29.70 29.61 17.39 17.29 1.20 57.82 1.00 2040 1750
## y0[22] 20.71 20.86 17.93 17.54 -9.24 49.87 1.00 1532 1971
## y0[23] 40.82 40.42 17.82 17.01 12.18 70.62 1.00 1796 1562
## y0[24] 45.77 46.12 18.00 17.66 16.56 76.43 1.00 1761 1836
## y0[25] 46.14 45.97 18.21 17.97 16.69 76.09 1.00 1846 1930
## y0[26] 43.59 43.90 17.77 17.39 14.22 72.23 1.00 1723 1871
## y0[27] 21.05 21.27 18.16 18.64 -8.34 50.24 1.00 1418 1892
## y0[28] 18.12 18.38 18.66 17.87 -12.75 49.07 1.00 1235 1857
## y0[29] 21.88 21.83 18.07 17.69 -7.16 51.98 1.00 1488 1966
## y0[30] 24.81 25.07 17.51 17.30 -5.84 53.77 1.00 1736 1736
## y0[31] 38.87 38.76 17.40 16.60 10.34 66.48 1.00 1770 1766
## y0[32] 39.88 39.49 17.59 17.99 11.82 68.76 1.00 2175 2015
## y0[33] 34.05 34.10 17.65 17.53 4.83 63.03 1.00 2083 1908
## y0[34] 22.93 22.66 17.93 17.80 -6.49 52.50 1.00 2074 1895
## y0[35] 41.96 42.07 17.96 17.70 13.19 71.52 1.00 2087 1970
## y0[36] 23.59 23.31 17.85 17.67 -5.25 53.40 1.00 1818 1998
## y0[37] 35.62 35.67 17.84 17.83 6.06 64.08 1.00 1961 2012
## y0[38] 43.46 42.78 17.98 17.73 15.21 74.84 1.00 2023 1674
## y0[39] 28.20 27.80 17.19 17.86 0.89 57.34 1.00 1863 2004
## y0[40] 30.84 30.69 17.71 17.80 2.75 59.08 1.00 1978 1930
## y0[41] 15.19 15.30 18.31 17.65 -15.93 45.09 1.00 1142 1642
## y0[42] 42.80 43.05 17.87 17.48 13.60 72.10 1.00 1981 1731
## y0[43] 22.94 22.43 17.76 17.52 -5.66 53.22 1.00 1785 2053
## y0[44] 25.70 25.73 17.56 16.95 -2.78 55.40 1.00 1719 1971
## y0[45] 37.84 38.35 18.08 17.44 8.63 67.26 1.00 1821 1787
## y0[46] 41.34 41.54 17.82 17.38 12.41 70.54 1.00 1994 1820
## y0[47] 40.67 41.20 17.29 16.97 12.87 69.87 1.00 1943 1709
## y0[48] 26.72 27.02 17.90 18.13 -2.21 56.44 1.00 1881 1878
## y0[49] 34.55 34.64 17.71 17.60 5.14 63.59 1.00 2144 2009
## y0[50] 22.29 22.64 17.46 17.94 -5.89 51.05 1.00 1804 1788
## m1[1] 14.88 14.91 5.75 5.70 4.93 23.57 1.02 216 256
## m1[2] 18.36 18.44 4.85 4.78 9.98 25.68 1.02 217 248
## m1[3] 21.84 21.90 4.01 3.94 14.96 27.88 1.02 223 245
## m1[4] 25.31 25.38 3.27 3.15 19.58 30.37 1.01 245 297
## m1[5] 28.79 28.78 2.70 2.63 24.21 33.01 1.01 318 487
## m1[6] 32.27 32.26 2.44 2.38 28.18 36.26 1.00 605 860
## m1[7] 35.74 35.70 2.58 2.48 31.41 39.98 1.00 1437 1347
## m1[8] 39.22 39.32 3.07 2.85 34.09 44.19 1.00 1243 1357
## m1[9] 42.70 42.82 3.76 3.57 36.35 48.79 1.01 798 1325
## m1[10] 46.17 46.30 4.58 4.40 38.39 53.70 1.01 578 1338
## y1[1] 14.72 14.61 18.42 18.70 -15.27 44.10 1.00 1537 1969
## y1[2] 18.77 19.14 18.18 17.70 -12.25 49.10 1.00 1229 1843
## y1[3] 22.28 21.90 18.27 17.93 -7.23 52.59 1.00 1724 1843
## y1[4] 25.33 25.71 17.97 17.62 -5.22 54.43 1.01 1256 1766
## y1[5] 29.05 28.70 17.41 17.04 0.81 57.79 1.00 1904 1988
## y1[6] 31.58 31.67 17.69 17.50 3.61 60.51 1.00 1963 1822
## y1[7] 35.99 36.43 17.54 16.97 7.30 64.84 1.00 2039 1856
## y1[8] 38.91 38.95 17.95 17.14 9.43 68.88 1.00 2050 2099
## y1[9] 42.30 42.42 17.99 18.15 12.97 72.76 1.00 1911 1974
## y1[10] 46.50 46.16 18.14 17.56 16.45 77.16 1.00 1821 1873
estimate as independent class
all b0l,b1l do not have a distribution
class c=1~k ex12-1.stan
yc~N(b0c+b1c*x,s)
data{
int N;
int K;
vector[N] x;
vector[N] y;
array[N] int c;
}
parameters{
vector[K] b0;
vector[K] b1;
real<lower=0> s;
}
model{
for(i in 1:N)
y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0[c[i]]+b1[c[i]]*x[i];
y0[i]=normal_rng(m0[i],s);
}
}mdl=cmdstan_model('./ex12-1.stan')
data=list(N=n,K=k,x=x,y=y,c=factor(c))
mle=mdl$optimize(data=data)## Initial log joint probability = -3454.57
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 -98.2382 0.0980673 0.370953 1 1 123
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 -98.2028 0.00186389 0.00452686 1 1 235
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 211 -98.2028 7.21993e-05 0.00177662 0.956 0.956 247
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -98.20
## b0[1] 5.71
## b0[2] -0.30
## b0[3] 15.80
## b0[4] 8.88
## b0[5] 12.25
## b0[6] 13.78
## b0[7] 16.90
## b0[8] 9.46
## b0[9] 8.57
## b1[1] 0.84
## b1[2] 4.24
## b1[3] -0.38
## b1[4] 4.18
## b1[5] 0.18
## b1[6] 2.01
## b1[7] 2.02
## b1[8] 2.61
## b1[9] 2.73
## s 4.32
## m0[1] 22.32
## m0[2] 35.32
## m0[3] 48.95
## m0[4] 29.75
## m0[5] 15.66
## m0[6] 15.25
## m0[7] 56.51
## m0[8] 10.16
## m0[9] 14.15
## m0[10] 48.87
## m0[11] 8.22
## m0[12] 14.37
## m0[13] 54.36
## m0[14] 37.45
## m0[15] 25.93
## m0[16] 26.58
## m0[17] 20.65
## m0[18] 14.26
## m0[19] 51.53
## m0[20] 42.89
## m0[21] 50.09
## m0[22] 21.82
## m0[23] 54.06
## m0[24] 53.20
## m0[25] 61.34
## m0[26] 15.56
## m0[27] 19.14
## m0[28] 21.31
## m0[29] 13.80
## m0[30] 26.67
## m0[31] 50.37
## m0[32] 46.06
## m0[33] 41.22
## m0[34] 24.19
## m0[35] 71.75
## m0[36] 13.31
## m0[37] 10.72
## m0[38] 84.94
## m0[39] 36.16
## m0[40] 53.51
## m0[41] 5.90
## m0[42] 9.14
## m0[43] 21.29
## m0[44] 28.54
## m0[45] 17.83
## m0[46] 19.32
## m0[47] 50.17
## m0[48] 32.03
## m0[49] 42.04
## m0[50] 26.93
## y0[1] 26.33
## y0[2] 36.61
## y0[3] 45.89
## y0[4] 28.22
## y0[5] 19.82
## y0[6] 14.93
## y0[7] 54.59
## y0[8] 11.15
## y0[9] 14.70
## y0[10] 53.16
## y0[11] 3.03
## y0[12] 19.77
## y0[13] 60.89
## y0[14] 41.39
## y0[15] 29.59
## y0[16] 23.59
## y0[17] 24.79
## y0[18] 13.43
## y0[19] 53.13
## y0[20] 48.38
## y0[21] 45.98
## y0[22] 16.53
## y0[23] 47.89
## y0[24] 52.90
## y0[25] 68.99
## y0[26] 13.91
## y0[27] 17.29
## y0[28] 18.32
## y0[29] 7.09
## y0[30] 24.85
## y0[31] 50.87
## y0[32] 44.12
## y0[33] 53.86
## y0[34] 23.39
## y0[35] 66.00
## y0[36] 13.44
## y0[37] 13.45
## y0[38] 84.43
## y0[39] 38.67
## y0[40] 60.75
## y0[41] 9.40
## y0[42] 3.88
## y0[43] 19.34
## y0[44] 29.15
## y0[45] 18.83
## y0[46] 12.34
## y0[47] 47.24
## y0[48] 29.06
## y0[49] 46.98
## y0[50] 39.90mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 1.6 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 1.6 seconds.
## Chain 3 finished in 1.6 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 1.6 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.6 seconds.
## Total execution time: 1.8 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -108.74 -108.29 3.82 3.69 -115.91 -103.30 1.00 570 818
## b0[1] 5.79 5.77 5.54 5.38 -3.22 14.84 1.01 2301 937
## b0[2] -0.36 -0.64 8.67 8.56 -14.09 14.14 1.00 1270 1262
## b0[3] 15.79 15.83 4.64 4.51 8.35 23.54 1.00 3237 1460
## b0[4] 9.04 8.55 11.58 11.35 -8.82 28.65 1.01 713 829
## b0[5] 12.46 12.52 7.68 7.97 -0.30 24.84 1.00 1328 1402
## b0[6] 13.84 13.90 5.56 5.37 4.75 23.10 1.00 3017 1255
## b0[7] 16.94 16.89 3.72 3.53 10.71 23.02 1.00 3335 1377
## b0[8] 9.79 9.65 7.31 7.14 -2.12 22.15 1.00 1840 1355
## b0[9] 8.57 8.53 3.94 3.90 2.15 15.33 1.00 2956 987
## b1[1] 0.83 0.84 0.39 0.39 0.17 1.45 1.00 2074 1446
## b1[2] 4.25 4.26 0.71 0.70 3.09 5.39 1.00 1335 1421
## b1[3] -0.39 -0.38 0.34 0.34 -0.95 0.17 1.00 2334 1594
## b1[4] 4.16 4.22 0.86 0.84 2.71 5.48 1.01 723 918
## b1[5] 0.18 0.17 0.63 0.65 -0.82 1.22 1.00 1360 1553
## b1[6] 2.01 2.00 0.42 0.40 1.32 2.71 1.00 2784 1626
## b1[7] 2.03 2.03 0.34 0.33 1.48 2.59 1.00 2861 1417
## b1[8] 2.59 2.59 0.66 0.65 1.49 3.68 1.00 1904 1333
## b1[9] 2.73 2.74 0.32 0.32 2.18 3.25 1.00 2676 1590
## s 5.61 5.57 0.71 0.68 4.57 6.88 1.00 1035 1337
## m0[1] 22.36 22.38 2.96 2.76 17.47 27.26 1.00 3401 1376
## m0[2] 35.42 35.37 2.55 2.49 31.36 39.61 1.00 1960 1740
## m0[3] 48.99 49.07 2.31 2.25 45.11 52.77 1.00 1915 1417
## m0[4] 29.80 29.81 2.13 2.02 26.23 33.27 1.00 3020 1521
## m0[5] 15.65 15.67 4.53 4.38 8.37 23.25 1.00 3249 1460
## m0[6] 15.26 15.18 3.27 3.20 9.79 20.75 1.00 2952 1165
## m0[7] 56.61 56.56 3.26 3.21 51.40 62.10 1.00 2003 1656
## m0[8] 10.12 10.07 2.58 2.62 5.98 14.33 1.00 1618 1149
## m0[9] 14.29 14.33 2.54 2.54 10.08 18.43 1.00 1740 1531
## m0[10] 48.95 48.96 2.77 2.73 44.33 53.44 1.01 2130 1365
## m0[11] 8.16 8.20 3.65 3.66 2.19 14.17 1.00 1585 1619
## m0[12] 14.50 14.54 2.44 2.43 10.49 18.51 1.00 1863 1505
## m0[13] 54.30 54.17 5.17 5.08 46.07 62.61 1.00 2046 1500
## m0[14] 37.54 37.55 2.54 2.47 33.51 41.74 1.00 1953 1698
## m0[15] 25.95 25.98 2.36 2.33 22.10 29.76 1.00 2841 1377
## m0[16] 26.75 26.72 3.57 3.49 20.89 32.87 1.00 1838 1443
## m0[17] 20.63 20.66 3.60 3.40 14.75 26.61 1.00 1967 1418
## m0[18] 14.39 14.42 2.46 2.47 10.31 18.47 1.00 1804 1439
## m0[19] 51.61 51.56 3.12 3.13 46.49 56.70 1.01 2186 1429
## m0[20] 42.96 42.98 2.13 2.08 39.42 46.40 1.01 2056 1167
## m0[21] 50.05 50.00 4.19 4.03 43.40 57.16 1.00 946 1146
## m0[22] 21.87 21.77 4.20 4.08 15.10 28.91 1.00 2983 1140
## m0[23] 54.11 54.16 2.70 2.54 49.54 58.56 1.01 1989 1338
## m0[24] 53.17 53.21 4.38 4.33 45.96 60.45 1.00 2113 1677
## m0[25] 61.39 61.43 3.36 3.22 55.76 66.93 1.01 2124 1551
## m0[26] 15.65 15.64 4.83 4.75 7.88 23.49 1.00 1643 1540
## m0[27] 19.15 19.09 2.91 2.87 14.30 23.95 1.00 2924 1260
## m0[28] 21.35 21.37 3.10 2.87 16.21 26.44 1.00 3423 1377
## m0[29] 13.78 13.81 3.22 3.16 8.60 19.03 1.00 3110 1356
## m0[30] 26.84 26.80 3.56 3.48 20.98 32.93 1.00 1839 1421
## m0[31] 50.41 50.47 2.41 2.30 46.32 54.33 1.00 1931 1379
## m0[32] 46.05 46.07 3.37 3.29 40.47 51.64 1.00 2046 1445
## m0[33] 41.29 41.33 2.00 1.96 38.03 44.53 1.01 2081 1294
## m0[34] 24.23 24.14 3.85 3.77 18.01 30.61 1.00 2943 1138
## m0[35] 71.90 71.99 4.94 4.98 63.60 80.01 1.00 1731 1487
## m0[36] 13.49 13.54 4.34 4.56 6.25 20.53 1.00 1356 1575
## m0[37] 10.67 10.68 2.42 2.47 6.80 14.64 1.00 1701 1366
## m0[38] 84.73 84.90 5.51 5.21 75.47 93.44 1.00 1106 1345
## m0[39] 36.20 36.21 3.56 3.48 30.27 42.20 1.00 1464 1422
## m0[40] 53.45 53.37 3.78 3.66 47.40 59.87 1.00 1075 1224
## m0[41] 5.98 5.97 5.46 5.36 -2.92 14.92 1.01 2304 937
## m0[42] 9.08 9.12 3.07 3.11 4.15 14.12 1.00 1576 1291
## m0[43] 21.31 21.27 2.73 2.70 16.78 25.79 1.00 2893 1259
## m0[44] 28.56 28.51 4.42 4.28 21.40 36.00 1.00 1341 1331
## m0[45] 17.84 17.83 2.94 2.85 12.98 22.74 1.00 2003 1554
## m0[46] 19.32 19.34 3.24 3.05 13.95 24.70 1.00 1970 1573
## m0[47] 50.24 50.21 2.94 2.91 45.40 55.02 1.01 2160 1338
## m0[48] 32.08 32.08 1.96 1.85 28.88 35.18 1.00 2826 1432
## m0[49] 42.07 42.05 1.96 1.86 38.92 45.20 1.00 1975 1293
## m0[50] 26.97 27.00 2.41 2.30 22.95 30.93 1.00 3145 1520
## y0[1] 22.38 22.34 6.18 6.04 12.18 32.72 1.00 2115 2055
## y0[2] 35.52 35.66 5.98 5.99 25.74 45.32 1.00 2192 1972
## y0[3] 48.90 49.10 6.06 5.99 38.75 58.79 1.00 2036 1996
## y0[4] 29.67 29.86 6.19 6.06 19.22 40.16 1.00 2252 1845
## y0[5] 15.46 15.48 7.31 7.09 3.06 27.56 1.00 2363 1685
## y0[6] 15.18 15.26 6.67 6.47 4.19 26.03 1.00 2347 2098
## y0[7] 56.52 56.51 6.52 6.53 45.52 67.42 1.00 1943 1892
## y0[8] 10.01 10.25 6.25 6.15 -0.11 20.33 1.00 1968 1916
## y0[9] 14.16 14.19 6.20 6.20 4.18 24.50 1.00 2091 1952
## y0[10] 49.08 49.06 6.46 6.45 38.57 59.71 1.00 2116 1598
## y0[11] 8.31 8.21 6.75 6.76 -2.58 19.51 1.00 1930 1787
## y0[12] 14.46 14.45 6.08 5.90 4.31 24.48 1.00 1898 1881
## y0[13] 54.28 54.35 7.82 7.59 41.47 66.65 1.00 1883 1772
## y0[14] 37.35 37.48 6.26 5.88 26.89 47.68 1.00 2057 2054
## y0[15] 25.70 25.62 6.02 5.79 15.92 35.56 1.00 1892 1726
## y0[16] 26.68 26.71 6.78 6.81 15.42 37.30 1.00 1877 1836
## y0[17] 20.49 20.31 6.88 6.88 9.46 31.83 1.00 2228 1748
## y0[18] 14.44 14.40 6.22 6.11 4.08 24.64 1.00 1902 1688
## y0[19] 51.50 51.69 6.47 6.41 40.97 62.02 1.00 1776 1634
## y0[20] 43.05 42.99 6.00 5.67 33.24 52.97 1.00 2035 1892
## y0[21] 50.21 50.29 7.02 6.99 38.80 61.73 1.00 1313 1582
## y0[22] 21.77 21.89 7.12 7.06 9.92 33.11 1.00 2380 1892
## y0[23] 54.15 54.35 6.24 6.21 43.32 63.75 1.00 1928 1788
## y0[24] 53.16 53.17 7.18 7.09 41.58 65.32 1.00 2069 1930
## y0[25] 61.57 61.58 6.54 6.30 50.90 72.28 1.00 2110 1894
## y0[26] 15.63 15.63 7.30 7.40 3.91 28.20 1.00 1900 1857
## y0[27] 19.06 19.01 6.25 5.94 8.88 29.03 1.00 2171 1999
## y0[28] 21.65 21.47 6.56 6.41 10.98 32.60 1.00 2220 1971
## y0[29] 13.76 13.54 6.51 6.38 3.53 24.65 1.00 1853 1791
## y0[30] 26.77 26.70 6.72 6.43 15.55 37.89 1.00 1978 1829
## y0[31] 50.45 50.46 6.25 6.06 40.09 61.17 1.00 2020 1923
## y0[32] 46.07 46.13 6.56 6.54 35.24 56.67 1.00 2031 1832
## y0[33] 41.26 41.35 6.08 5.95 31.36 51.02 1.00 2128 1931
## y0[34] 23.93 23.86 6.85 6.65 12.88 35.06 1.00 2397 1905
## y0[35] 71.93 71.92 7.64 7.50 59.56 84.53 1.00 1887 1880
## y0[36] 13.43 13.53 7.03 7.00 2.18 24.79 1.00 1770 1852
## y0[37] 10.59 10.47 6.12 6.15 0.60 20.92 1.00 1976 2014
## y0[38] 84.78 85.15 7.91 7.73 71.60 97.20 1.00 1487 1668
## y0[39] 36.30 36.25 6.72 6.44 25.36 47.29 1.00 1985 1931
## y0[40] 53.30 53.23 6.85 6.56 41.89 64.80 1.00 1512 1708
## y0[41] 5.99 6.15 7.93 7.94 -6.78 18.88 1.00 2044 1710
## y0[42] 9.25 9.28 6.38 6.26 -1.22 19.44 1.00 1570 1879
## y0[43] 21.42 21.22 6.18 6.17 11.81 31.56 1.00 2062 1974
## y0[44] 28.47 28.41 7.28 7.19 16.51 41.03 1.00 1809 1819
## y0[45] 17.92 17.92 6.30 6.18 7.54 28.40 1.00 1968 1799
## y0[46] 19.24 19.12 6.46 6.35 8.51 29.99 1.00 1904 1750
## y0[47] 50.53 50.43 6.31 6.29 40.28 61.10 1.00 2129 1449
## y0[48] 31.88 31.88 5.85 5.83 21.97 41.14 1.00 1936 1777
## y0[49] 41.97 41.91 5.89 5.79 32.24 51.91 1.00 2096 2057
## y0[50] 26.79 26.74 6.02 6.07 17.09 36.57 1.00 2015 2014
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,shape=c,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
#geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
#geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
#geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
#geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
qplot(x,y,facets=~c,size=I(2),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
estimate as class have relation
all b0l,b1l have a distribution
class c=1~k
b0c~N(b00,sb0)
b1c~N(b10,sb1)
yc~N(b0c+b1c*x,s)ex12-2.stan
class have relation
data{
int N;
int K;
vector[N] x;
vector[N] y;
array[N] int c;
}
parameters{
real b00;
real<lower=0,upper=100> sb0;
vector[K] b0;
real b10;
real<lower=0,upper=100> sb1;
vector[K] b1;
real<lower=0,upper=100> s;
}
model{
b0~normal(b00,sb0);
b1~normal(b10,sb1);
for(i in 1:N)
y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0[c[i]]+b1[c[i]]*x[i];
y0[i]=normal_rng(m0[i],s);
}
}mdl=cmdstan_model('./ex12-2.stan')
data=list(N=n,K=k,x=x,y=y,c=factor(c))
mle=mdl$optimize(data=data)## Initial log joint probability = -298.118
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 -29.1277 0.280031 1.12786e+06 0.6334 0.6334 166
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 163 0.0962256 3.18287e-05 1.38466e+06 0.8781 0.8781 256
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.try(print(mle))## variable estimate
## lp__ 0.10
## b00 0.92
## sb0 0.00
## b0[1] 0.92
## b0[2] 0.92
## b0[3] 0.92
## b0[4] 0.92
## b0[5] 0.92
## b0[6] 0.92
## b0[7] 0.92
## b0[8] 0.92
## b0[9] 0.92
## b10 2.14
## sb1 1.07
## b1[1] 1.88
## b1[2] 4.59
## b1[3] 1.92
## b1[4] 4.38
## b1[5] 2.30
## b1[6] 1.67
## b1[7] 2.94
## b1[8] 2.73
## b1[9] 2.23
## s 6.31
## m0[1] 8.81
## m0[2] 27.98
## m0[3] 33.94
## m0[4] 19.61
## m0[5] 1.62
## m0[6] 6.39
## m0[7] 62.47
## m0[8] 29.01
## m0[9] 24.78
## m0[10] 47.43
## m0[11] 38.71
## m0[12] 27.55
## m0[13] 47.92
## m0[14] 30.22
## m0[15] 15.12
## m0[16] 18.83
## m0[17] 34.60
## m0[18] 26.14
## m0[19] 51.29
## m0[20] 38.73
## m0[21] 44.10
## m0[22] 7.61
## m0[23] 38.12
## m0[24] 33.72
## m0[25] 44.06
## m0[26] 42.55
## m0[27] 9.56
## m0[28] 7.33
## m0[29] 10.90
## m0[30] 18.93
## m0[31] 35.10
## m0[32] 27.78
## m0[33] 36.30
## m0[34] 9.58
## m0[35] 78.98
## m0[36] 14.33
## m0[37] 26.26
## m0[38] 80.61
## m0[39] 40.42
## m0[40] 47.68
## m0[41] 1.35
## m0[42] 34.14
## m0[43] 11.33
## m0[44] 32.17
## m0[45] 28.25
## m0[46] 31.62
## m0[47] 49.31
## m0[48] 22.93
## m0[49] 28.28
## m0[50] 15.51
## y0[1] 13.94
## y0[2] 30.95
## y0[3] 36.00
## y0[4] 18.97
## y0[5] 6.09
## y0[6] 10.41
## y0[7] 63.29
## y0[8] 23.04
## y0[9] 21.54
## y0[10] 61.53
## y0[11] 36.39
## y0[12] 21.56
## y0[13] 46.39
## y0[14] 37.40
## y0[15] 16.80
## y0[16] 16.12
## y0[17] 24.91
## y0[18] 26.44
## y0[19] 56.45
## y0[20] 37.97
## y0[21] 46.34
## y0[22] 4.63
## y0[23] 38.75
## y0[24] 43.31
## y0[25] 44.43
## y0[26] 45.71
## y0[27] 15.61
## y0[28] 12.90
## y0[29] 9.94
## y0[30] 14.84
## y0[31] 22.51
## y0[32] 33.20
## y0[33] 37.90
## y0[34] 22.00
## y0[35] 80.42
## y0[36] 7.74
## y0[37] 19.19
## y0[38] 84.07
## y0[39] 33.08
## y0[40] 36.25
## y0[41] 10.50
## y0[42] 40.71
## y0[43] 14.35
## y0[44] 27.35
## y0[45] 27.76
## y0[46] 21.52
## y0[47] 38.06
## y0[48] 20.08
## y0[49] 27.12
## y0[50] 27.27mcmc=goMCMC(mdl,data,wrm=1000,smp=2000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 1 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 1 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 2 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 2 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 3 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 3 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 4 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 4 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 1 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 2 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 3 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 4 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 1 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 1 finished in 1.3 seconds.
## Chain 2 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 3 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 4 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 2 finished in 1.4 seconds.
## Chain 3 finished in 1.3 seconds.
## Chain 4 finished in 1.3 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.3 seconds.
## Total execution time: 1.4 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -127.67 -127.82 6.45 6.14 -138.27 -116.43 1.00 484 284
## b00 11.51 11.57 2.36 2.20 7.67 15.25 1.00 2577 3014
## sb0 3.58 3.04 2.47 2.13 0.75 8.24 1.00 525 273
## b0[1] 9.66 10.03 3.62 3.30 3.04 14.89 1.00 2747 4324
## b0[2] 10.11 10.48 4.07 3.47 2.75 15.98 1.00 3137 4413
## b0[3] 12.64 12.52 3.03 2.80 7.87 17.75 1.00 4199 4841
## b0[4] 12.06 11.97 4.24 3.39 5.59 19.13 1.00 4565 3786
## b0[5] 11.19 11.34 3.74 3.32 4.86 16.87 1.00 3977 4723
## b0[6] 12.29 12.19 3.19 2.94 7.30 17.56 1.00 4372 4316
## b0[7] 13.99 13.75 2.97 2.82 9.45 19.19 1.00 2132 2293
## b0[8] 11.29 11.44 3.49 3.06 5.32 16.82 1.00 4323 4076
## b0[9] 10.51 10.64 2.84 2.69 5.60 14.92 1.00 3084 3632
## b10 1.95 1.95 0.61 0.54 0.97 2.92 1.00 4096 2949
## sb1 1.70 1.59 0.57 0.45 1.04 2.74 1.00 2302 1311
## b1[1] 0.61 0.60 0.29 0.28 0.13 1.10 1.00 2699 940
## b1[2] 3.39 3.38 0.39 0.36 2.78 4.06 1.00 3045 1514
## b1[3] -0.16 -0.15 0.26 0.25 -0.58 0.27 1.00 2179 1097
## b1[4] 3.90 3.90 0.38 0.35 3.29 4.52 1.00 3926 3895
## b1[5] 0.30 0.29 0.36 0.33 -0.27 0.88 1.00 4740 4861
## b1[6] 2.11 2.11 0.29 0.28 1.63 2.58 1.00 4518 5289
## b1[7] 2.25 2.26 0.28 0.26 1.77 2.68 1.00 2342 2776
## b1[8] 2.44 2.43 0.37 0.34 1.83 3.04 1.00 4995 4857
## b1[9] 2.58 2.58 0.25 0.25 2.18 3.00 1.00 3815 4565
## s 5.47 5.40 0.67 0.64 4.50 6.70 1.00 2562 1449
## m0[1] 20.01 19.85 2.39 2.28 16.35 24.19 1.00 2343 3032
## m0[2] 35.42 35.44 2.25 2.18 31.71 39.15 1.00 5823 4865
## m0[3] 48.74 48.72 2.21 2.15 45.07 52.39 1.00 6963 4519
## m0[4] 28.26 28.19 1.80 1.73 25.44 31.32 1.00 3449 5194
## m0[5] 12.58 12.46 2.96 2.74 7.87 17.59 1.00 4227 4791
## m0[6] 16.84 16.95 2.38 2.23 12.75 20.48 1.00 3135 3448
## m0[7] 55.57 55.58 3.04 3.04 50.54 60.57 1.00 5032 1637
## m0[8] 10.32 10.33 2.54 2.46 6.18 14.44 1.00 2105 966
## m0[9] 14.27 14.26 2.22 2.15 10.63 17.89 1.00 5654 4662
## m0[10] 49.50 49.52 2.53 2.50 45.37 53.60 1.00 5780 5241
## m0[11] 9.51 9.51 3.49 3.34 3.89 15.25 1.00 1879 952
## m0[12] 14.62 14.61 2.37 2.27 10.75 18.52 1.00 6008 4650
## m0[13] 53.20 53.23 4.11 4.02 46.54 60.01 1.00 5880 5218
## m0[14] 37.42 37.44 2.38 2.28 33.50 41.32 1.00 5962 5026
## m0[15] 26.95 27.01 1.82 1.77 23.87 29.80 1.00 3694 3905
## m0[16] 27.26 27.29 2.12 2.04 23.78 30.65 1.00 4930 4995
## m0[17] 20.55 20.50 3.48 3.46 14.94 26.25 1.00 3457 3332
## m0[18] 14.44 14.43 2.28 2.21 10.69 18.19 1.00 5866 4765
## m0[19] 52.45 52.49 2.82 2.76 47.82 57.01 1.00 5016 5446
## m0[20] 42.86 42.86 1.99 1.90 39.60 46.08 1.00 8074 5164
## m0[21] 50.50 50.55 2.68 2.62 46.01 54.84 1.00 2976 1030
## m0[22] 20.72 20.64 2.51 2.27 16.76 24.95 1.00 4860 4403
## m0[23] 53.59 53.55 2.54 2.50 49.38 57.76 1.00 6697 4719
## m0[24] 53.63 53.65 4.18 4.04 46.73 60.47 1.00 5691 5324
## m0[25] 60.47 60.47 3.08 2.99 55.38 65.47 1.00 6116 6405
## m0[26] 16.56 16.46 3.98 3.88 10.06 23.18 1.00 5725 5864
## m0[27] 20.51 20.61 2.14 2.03 16.81 23.77 1.00 3235 3792
## m0[28] 18.89 18.70 2.50 2.37 15.06 23.22 1.00 2285 3081
## m0[29] 11.82 11.73 2.22 2.11 8.27 15.59 1.00 4459 4878
## m0[30] 27.35 27.38 2.12 2.05 23.86 30.74 1.00 4937 5011
## m0[31] 50.09 50.06 2.30 2.24 46.28 53.89 1.00 6901 4690
## m0[32] 46.15 46.12 3.36 3.21 40.54 51.68 1.00 5979 5008
## m0[33] 41.00 41.00 1.87 1.78 37.93 44.04 1.00 8174 5268
## m0[34] 23.20 23.13 2.39 2.20 19.41 27.24 1.00 5055 4138
## m0[35] 67.77 67.85 4.08 4.01 60.90 74.43 1.00 4162 1627
## m0[36] 12.92 12.96 2.36 2.17 8.93 16.60 1.00 4153 5037
## m0[37] 10.54 10.56 2.34 2.26 6.75 14.36 1.00 2252 967
## m0[38] 83.00 83.10 4.51 4.38 75.20 90.35 1.00 2951 976
## m0[39] 39.29 39.30 2.36 2.30 35.43 43.09 1.00 5058 5867
## m0[40] 53.68 53.75 2.75 2.69 49.07 58.16 1.00 2798 1018
## m0[41] 9.80 10.16 3.57 3.27 3.26 14.96 1.00 2758 4443
## m0[42] 9.89 9.86 3.01 2.92 5.00 14.76 1.00 1933 988
## m0[43] 22.56 22.65 2.02 1.93 19.08 25.67 1.00 3335 3958
## m0[44] 33.19 33.25 2.45 2.30 29.05 36.98 1.00 4399 5433
## m0[45] 18.50 18.45 2.86 2.81 13.90 23.30 1.00 3973 4705
## m0[46] 19.59 19.53 3.17 3.12 14.48 24.89 1.00 3689 4840
## m0[47] 50.94 50.97 2.67 2.63 46.55 55.26 1.00 5367 5460
## m0[48] 30.79 30.74 1.71 1.64 28.09 33.68 1.00 4333 5527
## m0[49] 42.20 42.20 1.86 1.81 39.13 45.27 1.00 6548 4363
## m0[50] 25.13 25.02 1.99 1.92 22.05 28.52 1.00 2840 4361
## y0[1] 19.94 19.87 6.05 5.85 10.00 29.94 1.00 6955 7166
## y0[2] 35.36 35.37 5.91 5.78 25.65 45.00 1.00 7888 7091
## y0[3] 48.69 48.74 6.00 5.95 38.90 58.53 1.00 8348 7735
## y0[4] 28.16 28.20 5.76 5.73 18.66 37.58 1.00 7904 7642
## y0[5] 12.65 12.70 6.34 6.26 2.21 22.90 1.00 6808 7459
## y0[6] 16.81 16.88 6.02 5.90 7.04 26.63 1.00 7310 7592
## y0[7] 55.52 55.51 6.24 6.04 45.16 65.61 1.00 8259 7526
## y0[8] 10.31 10.35 6.07 5.96 0.30 20.27 1.00 5976 7655
## y0[9] 14.33 14.36 5.88 5.79 4.61 23.92 1.00 8158 7779
## y0[10] 49.55 49.61 6.05 5.95 39.67 59.39 1.00 7730 7677
## y0[11] 9.44 9.46 6.49 6.42 -1.01 19.91 1.00 4661 4641
## y0[12] 14.57 14.53 5.98 5.84 4.94 24.52 1.00 7647 7719
## y0[13] 53.24 53.32 6.93 6.68 41.78 64.41 1.00 7089 7608
## y0[14] 37.44 37.45 5.96 5.86 27.78 47.22 1.00 7957 7632
## y0[15] 27.00 26.95 5.84 5.85 17.61 36.47 1.00 7640 7800
## y0[16] 27.27 27.25 5.82 5.80 17.68 36.79 1.00 7342 7930
## y0[17] 20.62 20.62 6.52 6.25 9.84 31.38 1.00 7453 7719
## y0[18] 14.36 14.45 6.01 5.84 4.24 24.25 1.00 8288 8057
## y0[19] 52.47 52.39 6.13 5.96 42.68 62.60 1.00 7836 7847
## y0[20] 42.88 42.86 5.84 5.75 33.25 52.41 1.00 8276 7855
## y0[21] 50.44 50.39 6.14 5.98 40.32 60.51 1.00 7520 7321
## y0[22] 20.79 20.72 6.01 5.92 11.10 30.75 1.00 7117 7280
## y0[23] 53.57 53.56 6.00 5.93 43.74 63.24 1.00 7770 6694
## y0[24] 53.62 53.52 6.87 6.70 42.32 64.88 1.00 7546 6710
## y0[25] 60.51 60.48 6.30 6.11 50.26 70.76 1.00 8320 7630
## y0[26] 16.59 16.56 6.81 6.63 5.52 28.03 1.00 7209 7013
## y0[27] 20.42 20.46 5.88 5.77 10.69 30.12 1.00 6954 7102
## y0[28] 18.96 18.99 6.05 6.08 9.06 28.82 1.00 6721 7176
## y0[29] 11.80 11.83 5.92 5.94 2.21 21.49 1.00 8079 7972
## y0[30] 27.27 27.28 5.89 5.73 17.52 36.94 1.00 7868 7731
## y0[31] 50.16 50.14 5.96 5.82 40.48 60.14 1.00 7698 7812
## y0[32] 46.14 46.22 6.41 6.37 35.83 56.73 1.00 7767 7435
## y0[33] 40.94 40.98 5.80 5.63 31.35 50.45 1.00 7887 7649
## y0[34] 23.12 23.10 6.06 5.93 13.34 33.13 1.00 7822 7161
## y0[35] 67.74 67.83 6.82 6.58 56.45 78.81 1.00 5690 6428
## y0[36] 13.00 13.01 5.98 5.84 3.10 22.81 1.00 6983 7602
## y0[37] 10.60 10.65 6.02 5.84 0.66 20.48 1.00 6856 7218
## y0[38] 83.02 83.07 7.16 7.07 71.14 94.68 1.00 4983 5998
## y0[39] 39.34 39.39 5.93 5.94 29.79 49.02 1.00 7688 7518
## y0[40] 53.63 53.68 6.19 6.14 43.34 63.67 1.00 6851 7788
## y0[41] 9.81 9.72 6.57 6.50 -0.93 20.56 1.00 4819 6683
## y0[42] 9.83 9.81 6.30 6.13 -0.45 20.23 1.00 5210 6413
## y0[43] 22.58 22.57 5.83 5.75 13.10 32.15 1.00 6605 7333
## y0[44] 33.30 33.31 6.09 6.03 23.29 43.33 1.00 7359 7834
## y0[45] 18.43 18.54 6.25 6.19 7.97 28.50 1.00 7375 6806
## y0[46] 19.56 19.59 6.34 6.26 9.13 30.03 1.00 6599 7050
## y0[47] 50.95 50.97 6.15 6.14 41.09 61.11 1.00 7560 7443
## y0[48] 30.83 30.88 5.74 5.75 21.33 40.17 1.00 7571 7651
## y0[49] 42.18 42.22 5.76 5.67 32.74 51.87 1.00 7699 7161
## y0[50] 25.16 25.14 5.86 5.84 15.49 34.87 1.00 6799 7173
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,shape=c,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
#geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
#geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
#geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
#geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
qplot(x,y,facets=~c,size=I(2),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
generalized linear mixed model
(X,y)i=1-n
b[b0,b1,...]
linear model y~N(Xb,s)
generalized linear model y~dist.(m=link(Xb),s)
fixed effect b0, b1,...
individual random effect b0+r0i r0~N(0,sr0), b1+r1i r1~N(0,sr1),...
class c=1-k
class effect b0+r0c r0~N(0,sr0), b1+r1c r1~N(0,sr1),...
for y=dist.(m,s)
random intercept model m=b0+r0i+b1*x, m=b0+r0c+b1*x
random coefficient model m=b0+(b1+r1i)*x, m=b0+(b1+r1c)*x
mixed model m=b0+r0i+(b1+r1i)*x, m=b0+r0c+(b1+r1c)*x
note
@ yi=b0+b1*xi+r0i is not useful, ri is included in s
@ yi=b0+b1*xi+r0c is useful, repeated observation can be treated by class effect
@ when appling Poisson dist.(y is larger than 0, error is larger at large x),
but variance is larger than mean (over dispersion),
random intercept model is usefuln=20
x=runif(n,0,10)
r0=rnorm(n,0,1)
y=rpois(n,exp(2+0.1*x+r0))
qplot(x,y)
ex13-0.stan
generalized linear model, poisson regression
data{
int N;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
}
model{
y~poisson_log(b0+b1*x);
}
generated quantities{
vector[N] m0;
array[N] int y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=poisson_log_rng(m0[i]);
}
}mdl=cmdstan_model('./ex13-0.stan')
data=list(N=n,x=x,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -3.4008e+06
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 18 495.512 0.00132099 0.0675192 1 1 48
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 495.51
## b0 2.23
## b1 0.09
## m0[1] 2.24
## m0[2] 2.25
## m0[3] 2.44
## m0[4] 2.30
## m0[5] 3.01
## m0[6] 2.82
## m0[7] 2.52
## m0[8] 2.94
## m0[9] 3.05
## m0[10] 2.76
## m0[11] 2.62
## m0[12] 2.52
## m0[13] 3.02
## m0[14] 2.95
## m0[15] 2.66
## m0[16] 2.55
## m0[17] 2.38
## m0[18] 2.54
## m0[19] 2.86
## m0[20] 2.35
## y0[1] 9.00
## y0[2] 8.00
## y0[3] 11.00
## y0[4] 14.00
## y0[5] 19.00
## y0[6] 17.00
## y0[7] 6.00
## y0[8] 11.00
## y0[9] 20.00
## y0[10] 18.00
## y0[11] 17.00
## y0[12] 13.00
## y0[13] 18.00
## y0[14] 20.00
## y0[15] 12.00
## y0[16] 18.00
## y0[17] 10.00
## y0[18] 17.00
## y0[19] 13.00
## y0[20] 14.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 494.48 494.83 1.02 0.73 492.41 495.46 1.00 639 984
## b0 2.23 2.23 0.12 0.12 2.02 2.43 1.01 380 404
## b1 0.09 0.09 0.02 0.02 0.05 0.12 1.00 435 655
## m0[1] 2.24 2.24 0.12 0.12 2.03 2.44 1.01 381 404
## m0[2] 2.25 2.25 0.12 0.12 2.05 2.45 1.01 381 404
## m0[3] 2.44 2.44 0.08 0.08 2.30 2.58 1.01 416 495
## m0[4] 2.30 2.30 0.11 0.11 2.12 2.48 1.01 384 417
## m0[5] 3.00 3.00 0.09 0.09 2.85 3.16 1.00 1129 1200
## m0[6] 2.82 2.82 0.07 0.06 2.72 2.93 1.00 1767 1477
## m0[7] 2.51 2.51 0.07 0.07 2.39 2.63 1.01 463 546
## m0[8] 2.94 2.94 0.08 0.08 2.80 3.08 1.00 1426 1332
## m0[9] 3.05 3.05 0.10 0.10 2.88 3.22 1.00 1009 1070
## m0[10] 2.76 2.76 0.06 0.06 2.66 2.86 1.00 1536 1302
## m0[11] 2.62 2.62 0.06 0.06 2.51 2.72 1.00 654 897
## m0[12] 2.51 2.51 0.07 0.07 2.39 2.63 1.01 463 546
## m0[13] 3.01 3.01 0.09 0.09 2.86 3.17 1.00 1104 1200
## m0[14] 2.95 2.95 0.08 0.08 2.81 3.09 1.00 1371 1351
## m0[15] 2.66 2.65 0.06 0.06 2.55 2.76 1.00 809 987
## m0[16] 2.55 2.55 0.07 0.07 2.44 2.67 1.00 511 720
## m0[17] 2.38 2.38 0.09 0.09 2.22 2.54 1.01 396 433
## m0[18] 2.54 2.54 0.07 0.07 2.42 2.66 1.00 493 650
## m0[19] 2.85 2.85 0.07 0.07 2.74 2.97 1.00 1742 1441
## m0[20] 2.35 2.35 0.10 0.10 2.18 2.51 1.01 390 423
## y0[1] 9.54 9.00 3.32 2.97 5.00 15.00 1.00 1279 1882
## y0[2] 9.56 9.00 3.31 2.97 5.00 16.00 1.00 1506 1568
## y0[3] 11.47 11.00 3.52 2.97 6.00 18.00 1.00 1478 1928
## y0[4] 9.99 10.00 3.29 2.97 5.00 16.00 1.00 1655 1773
## y0[5] 20.28 20.00 4.86 4.45 13.00 29.00 1.00 1644 1847
## y0[6] 16.72 17.00 4.23 4.45 10.00 24.00 1.00 1738 1880
## y0[7] 12.47 12.00 3.74 4.45 7.00 19.00 1.00 1741 1761
## y0[8] 18.94 19.00 4.66 4.45 12.00 27.00 1.00 1852 1945
## y0[9] 21.03 21.00 5.01 4.45 13.00 30.00 1.00 1803 1839
## y0[10] 15.76 16.00 4.12 4.45 9.00 23.00 1.00 1776 1845
## y0[11] 13.85 14.00 3.77 2.97 8.00 20.00 1.00 1970 2055
## y0[12] 12.46 12.00 3.68 2.97 7.00 19.00 1.00 1782 1982
## y0[13] 20.61 20.00 4.97 4.45 13.00 29.00 1.00 1761 1797
## y0[14] 19.10 19.00 4.70 4.45 12.00 27.00 1.00 2200 1998
## y0[15] 14.14 14.00 3.88 4.45 8.00 21.00 1.00 1957 1752
## y0[16] 12.80 13.00 3.64 2.97 7.00 19.00 1.00 1796 2006
## y0[17] 10.89 11.00 3.49 2.97 5.00 17.00 1.00 1427 1653
## y0[18] 12.52 12.00 3.73 2.97 7.00 19.00 1.00 1704 1960
## y0[19] 17.50 17.00 4.43 4.45 11.00 25.00 1.00 2163 1923
## y0[20] 10.57 10.00 3.25 2.97 5.00 16.00 1.00 1385 1821
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
ex13-1.stan
generalized linear mixed model, poisson regression
data{
int N;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
real<lower=0> sr0;
vector[N] r0;
}
model{
r0~normal(0,sr0);
for(i in 1:N)
y~poisson_log(b0+r0[i]+b1*x);
}
generated quantities{
vector[N] m0;
array[N] int y0;
for(i in 1:N){
m0[i]=b0+r0[i]+b1*x[i];
y0[i]=poisson_log_rng(m0[i]);
}
}mdl=cmdstan_model('./ex13-1.stan')
data=list(N=n,x=x,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -30405.3
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 10069.1 0.0268383 19129.8 0.2113 0.02113 145
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 10309.9 0.0331786 8.55093e+13 0.3351 0.3351 342
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 206 10310 0.00319902 1.03472e+12 0.681 0.681 349
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 10310.00
## b0 2.84
## b1 0.00
## sr0 0.00
## r0[1] 0.00
## r0[2] 0.00
## r0[3] 0.00
## r0[4] 0.00
## r0[5] 0.00
## r0[6] 0.00
## r0[7] 0.00
## r0[8] 0.00
## r0[9] 0.00
## r0[10] 0.00
## r0[11] 0.00
## r0[12] 0.00
## r0[13] 0.00
## r0[14] 0.00
## r0[15] 0.00
## r0[16] 0.00
## r0[17] 0.00
## r0[18] 0.00
## r0[19] 0.00
## r0[20] 0.00
## m0[1] 2.84
## m0[2] 2.84
## m0[3] 2.84
## m0[4] 2.84
## m0[5] 2.84
## m0[6] 2.84
## m0[7] 2.84
## m0[8] 2.84
## m0[9] 2.84
## m0[10] 2.84
## m0[11] 2.84
## m0[12] 2.84
## m0[13] 2.84
## m0[14] 2.84
## m0[15] 2.84
## m0[16] 2.84
## m0[17] 2.84
## m0[18] 2.84
## m0[19] 2.84
## m0[20] 2.84
## y0[1] 12.00
## y0[2] 15.00
## y0[3] 21.00
## y0[4] 12.00
## y0[5] 15.00
## y0[6] 26.00
## y0[7] 17.00
## y0[8] 14.00
## y0[9] 9.00
## y0[10] 17.00
## y0[11] 14.00
## y0[12] 19.00
## y0[13] 21.00
## y0[14] 15.00
## y0[15] 14.00
## y0[16] 21.00
## y0[17] 22.00
## y0[18] 17.00
## y0[19] 13.00
## y0[20] 9.00mcmc=goMCMC(mdl,data,wrm=500)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1 finished in 0.8 seconds.
## Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 3 finished in 0.8 seconds.
## Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 2 finished in 0.9 seconds.
## Chain 4 finished in 0.8 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.8 seconds.
## Total execution time: 1.0 seconds.seeMCMC(mcmc,'[m,r]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 9984.76 9983.98 13.47 13.68 9963.32 10008.50 1.02 76 92
## b0 2.23 2.23 0.03 0.03 2.19 2.28 1.00 702 722
## b1 0.09 0.09 0.00 0.00 0.08 0.09 1.00 784 743
## sr0 0.01 0.01 0.01 0.01 0.00 0.03 1.02 76 91
## r0[1] 0.00 0.00 0.02 0.01 -0.03 0.03 1.01 4126 787
## r0[2] 0.00 0.00 0.01 0.01 -0.03 0.02 1.01 2956 887
## r0[3] 0.00 0.00 0.02 0.01 -0.03 0.03 1.01 2947 724
## r0[4] 0.00 0.00 0.01 0.01 -0.03 0.02 1.00 3740 873
## r0[5] 0.00 0.00 0.02 0.01 -0.02 0.03 1.01 2896 714
## r0[6] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 3065 713
## r0[7] 0.00 0.00 0.02 0.01 -0.03 0.03 1.01 2504 735
## r0[8] 0.00 0.00 0.02 0.01 -0.02 0.02 1.00 3240 907
## r0[9] 0.00 0.00 0.01 0.01 -0.02 0.02 1.00 2565 955
## r0[10] 0.00 0.00 0.02 0.01 -0.02 0.02 1.00 2608 743
## r0[11] 0.00 0.00 0.02 0.01 -0.02 0.03 1.01 3916 828
## r0[12] 0.00 0.00 0.02 0.01 -0.02 0.02 1.01 3169 605
## r0[13] 0.00 0.00 0.02 0.01 -0.03 0.03 1.01 3166 778
## r0[14] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 3116 965
## r0[15] 0.00 0.00 0.02 0.01 -0.02 0.02 1.00 2852 909
## r0[16] 0.00 0.00 0.01 0.01 -0.02 0.02 1.00 3048 774
## r0[17] 0.00 0.00 0.01 0.01 -0.02 0.02 1.00 3239 1041
## r0[18] 0.00 0.00 0.02 0.01 -0.03 0.03 1.00 2791 694
## r0[19] 0.00 0.00 0.01 0.01 -0.03 0.02 1.01 3670 862
## r0[20] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 3273 959
## m0[1] 2.24 2.24 0.03 0.03 2.19 2.29 1.00 842 852
## m0[2] 2.25 2.25 0.03 0.03 2.20 2.30 1.00 844 998
## m0[3] 2.45 2.45 0.02 0.02 2.41 2.48 1.00 986 943
## m0[4] 2.30 2.30 0.03 0.03 2.25 2.35 1.01 881 1053
## m0[5] 3.01 3.01 0.02 0.03 2.97 3.05 1.00 1774 1533
## m0[6] 2.82 2.82 0.02 0.02 2.79 2.86 1.00 2348 1663
## m0[7] 2.52 2.52 0.02 0.02 2.48 2.55 1.00 1198 968
## m0[8] 2.94 2.94 0.02 0.02 2.90 2.98 1.00 2075 1350
## m0[9] 3.05 3.05 0.03 0.03 3.00 3.09 1.00 1509 1359
## m0[10] 2.76 2.76 0.02 0.02 2.73 2.79 1.00 2515 986
## m0[11] 2.62 2.62 0.02 0.02 2.59 2.65 1.01 1927 901
## m0[12] 2.52 2.52 0.02 0.02 2.48 2.55 1.00 1137 1189
## m0[13] 3.02 3.02 0.03 0.02 2.97 3.06 1.00 1898 1388
## m0[14] 2.95 2.95 0.02 0.02 2.91 2.99 1.00 1773 1269
## m0[15] 2.66 2.66 0.02 0.02 2.63 2.69 1.00 1972 1231
## m0[16] 2.56 2.56 0.02 0.02 2.52 2.59 1.00 1342 1128
## m0[17] 2.39 2.38 0.03 0.02 2.35 2.43 1.00 922 1036
## m0[18] 2.54 2.54 0.02 0.02 2.51 2.58 1.00 1111 1130
## m0[19] 2.86 2.86 0.02 0.02 2.82 2.89 1.00 2567 1308
## m0[20] 2.35 2.35 0.03 0.03 2.31 2.39 1.00 961 987
## y0[1] 9.44 9.00 2.99 2.97 5.00 15.00 1.00 2104 2052
## y0[2] 9.53 9.00 3.04 2.97 5.00 15.00 1.00 1947 2018
## y0[3] 11.50 11.00 3.40 2.97 6.00 17.00 1.00 1989 1787
## y0[4] 9.98 10.00 3.22 2.97 5.00 16.00 1.00 1898 2005
## y0[5] 20.23 20.00 4.62 4.45 13.00 28.00 1.00 2159 1847
## y0[6] 16.87 17.00 4.18 4.45 10.00 24.00 1.00 1887 1877
## y0[7] 12.36 12.00 3.65 2.97 7.00 19.00 1.00 2082 2040
## y0[8] 18.91 19.00 4.38 4.45 12.00 26.00 1.00 2031 2028
## y0[9] 21.25 21.00 4.64 4.45 14.00 29.00 1.00 1919 1885
## y0[10] 15.76 16.00 3.93 4.45 10.00 22.00 1.00 2012 2037
## y0[11] 13.79 14.00 3.71 4.45 8.00 20.00 1.00 2045 2075
## y0[12] 12.24 12.00 3.48 2.97 7.00 18.00 1.00 2077 2035
## y0[13] 20.44 20.00 4.56 4.45 13.00 28.00 1.00 1987 1818
## y0[14] 19.15 19.00 4.41 4.45 12.00 27.00 1.00 1997 2056
## y0[15] 14.21 14.00 3.73 2.97 8.00 21.00 1.00 2104 1893
## y0[16] 12.85 13.00 3.57 2.97 7.00 19.00 1.00 1981 1923
## y0[17] 10.80 11.00 3.34 2.97 6.00 17.00 1.00 1989 1937
## y0[18] 12.69 13.00 3.62 4.45 7.00 19.00 1.00 2127 1968
## y0[19] 17.40 17.00 4.21 4.45 11.00 25.00 1.00 2019 1915
## y0[20] 10.43 10.00 3.32 2.97 5.00 16.00 1.00 1962 1853
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')